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ABSTRACT 

Extragalactic jets are visualized as dynamic erruptive events modelled by time- 
dependent magnetohydrodynamic (MHD) equations. The jet structure comes 
through the temporally self-similar solutions in two-dimensional axisymmetric 
spherical geometry. The two-dimensional magnetic field is solved in the finite plasma 
pressure regime, or finite j3 regime, and it is described by an equation where plasma 
pressure plays the role of an eigenvalue. This allows a structure of magnetic lobes 
in space, among which the polar axis lobe is strongly peaked in intensity and col- 
limated in angular spread comparing to the others. For this reason, the polar lobe 
overwhelmes the other lobes, and a jet structure arises in the polar direction natu- 
rally. Furthermore, within each magnetic lobe in space, there are small secondary 
regions with closed two-dimensional field lines embedded along this primary lobe. 
In these embedded magnetic toroids, plasma pressure and mass density are much 
higher accordingly. These are termed as secondary plasmoids. The magnetic field 
lines in these secondary plasmoids circle in alternating sequence such that adjacent 
plasmoids have opposite field lines. In particular, along the polar primary lobe, 
such periodic plasmoid structure happens to be compatible with radio observations 
where islands of high radio intensities are mapped. 



Subject headings: Jets, Accretion Disks, Self-Similar MHD 
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1. Introduction 

CoUimated jets with high terminal velocities appear to be universal phenomena in astro- 
physics. They are often associated with young stellar objects, compact galactic objects, and 
active galactic nuclei [Livio 1997]. These jets are always accompanied by accretion disks on 
the equatorial plane. The magnetosphere of the accretion disk contains a plasma that follows 
the accretion of the materials in the disk. The plasma density and pressure get higher as the 
accretion approaches to the center which drive the magnetic field. The dynamics of this system 
was first described in magnetohydrodynamic (MHD) model by Blandford and Payne [1982]. In 
this landmark paper, jets are considered as a spatial structure in stationary state. Ideal MHD 
equations in cylindrical coordinates (r, 0, z) are solved for time independent solutions. The 
central mass M is replaced by a linear mass along the cylindrical axis. Self-similar solutions 
in space with a scale invariance z/r are sought. Such self-similar solutions are compatible to 
a Keplerian disk plasma rotation velocity field superimposed by an Alfvenic plasma velocity. 
The complete MHD instability spectrum with such Keplerian profile are analysed by Keppens, 
Casse, and Goedbloed [2002]. This model predicts that the magnetospheric disk plasma would 
be ejected towards the polar direction magnetocentrifugally should the magnetic field lines be 
at an angle less than 7r/3 or more that 27r/3 on the rz plane with respect to the outward radius 
of the disk. CoUimating action of this plasma outflow would be provided by the hoop force of 
the azimuthal magnetic field and its associated parallel current in a force-free configuration. 
The interaction of the magnetic field with the plasma disk generates a MHD Poynting flux [Fer- 
reira and Pelletier 1995] that can be converted into kinetic energy of the jet plasma [Zanni et. 
al. 2004]. Variants of jet formation model are proposed by Contopoulos and Lovelace [1994], 
Contopoulos [1995], and Cao and Spruit [1994]. This accretion-ejection model provides the 
basic framework of current investigations of accretion disk and jets as is reviewed by Balbus 
and Hawley [1998]. Dissipative MHD effects are examined by Casse and Ferreira [2000 a,b] 
and Casse [2004] , and relativistic jets are analyzed by Vlahakis [2004] . These stationary state 
analytic studies are often complemented by numerical works in cylindrical geometry to simulate 
time evolutions of the jets [Ustyugova et. al. 1995, Ouyed and Pudritz 1997, Krasnopolski et. 
al. 1999], and the disk-jet system [Matsumoto et. al. 1996, Kato et. al. 2002]. 

Here, instead of a stationary model, we take the view that jets are a time dependent spatial 
structure. What we see is only a snapshot of their state at this particular moment. Due to 
their galactic dimensions, the time scale of these structures is believed to be extraordinarily 
large which gives the impression of a stationary structure. This implies that jets are results 
of an erruption originating from the galactic nucleus. They could be dissipated in time before 
another erruption takes place due to pressure built-up from accretion. Or one erruption could 
be superimposed on an earlier event. To model the jet system, we will do a self-similar analysis 
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on the full time-dependent ideal MHD equations in spherical coordinates (r, 9^ 0) with a mass 
Af at the nucleus. In particular, we consider axisymmetric solutions. This type of self-similar 
solutions were pioneered by Low [1982a,b, 1984a,b] for astrophysical applications and solar 
corona mass ejections with pure radial plasma velocity flow. Variants of these solutions include 
cases where the plasma domain lies outside the mass M such as interplanetary magnetic ropes in 
one-dimensional [Osherovich, Farrugia, and Burlaga 1993, 1995] and in two-dimensional [Tsui 
and Tavares 2005] cylindrical geometry, interplanetary magnetic clouds [Tsui 2006a], and also 
atmospheric ball lightnings [Tsui 2006b] in spherical geometry. In these descriptions where M 
is outside the spherical domain of interest, the magnetic field is axisymmetric force-free and 
contains regions of closed field lines, while the plasma is spherically symmetric decoupled from 
the magnetic field. 

For the present case of extragalactic jets described by the mechanism of accretion-ejection, 
we will follow the self-similar solutions of Low with the polytropic index 7 = 4/3, but with 
particular emphasis on the finite plasma pressure. This current approach differs from the 
Keplerian disk plasma in that the radial flow is not tied to the Alfven velocity as a priori. In 
this dynamic model, we consider jets as a manifestation of mass ejection on a galactic scale. The 
plasma pressure, in this self-similar MHD model, proves to have an important role in coUimating 
the magnetic fields and the jet plasmas. The time evolution function gives a dynamic description 
of the high radial flow especially in the jets. The self-similar solutions, that converge at the 
center and at inflnity, give small regions of closed axisymmetric two-dimensional magnetic fleld 
lines where plasma density and pressure are much higher. These regions along the jets could 
correspond to the high intensity islands in radio frequency maps. 
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2. Self-Similar MHD 

Following the accretion-ejection classical model, we also use the MHD equations to de- 
scribe the plasma. Nevertheless, we retain the time dependence to write 



^ + V-M = 0, (1) 

+ (v-V)v} ^ JxB-Vp- p^r , (2) 

— ^ -VxE ^Vx{vxB) , (3) 

VxB = i^J , (4) 

V ■ 5 = , (5) 

^{-) + {v-V){^) = . (6) 



Here, p is the mass density, v is the bulk velocity, J is the current density, B is the magnetic 
field, "p is the plasma pressure, p is the free space permeability, and 7 is the polytropic index. 
For the bulk velocity, it is consisted of a radial and an azimuthal component to model the 
plasma outflow and the disk rotation. To be compatible to the physical situation, the meridian 
velocity is taken to be null. For the radial component ff, we seek self-similar solutions where 
the time evolution is described by the dimensionless evolution function yii). The radial profile 
is time invariant in terms of the radial label ri = r{t)/y{t), which has the dimension of r, such 
that 1] is independent of time, and corresponds to the Lagrangian radial position attached to 
a fixed fiuid element. The label rj corresponds to the radial positions of the initial self-similar 
configuration that expands in time with radial Eulerian positions r{t) = f]y{t). The velocity 
can then be written as 

^ = {^5' 0,^4 . (7) 



We consider a two-dimensional case with azimuthal symmetry in 0. In this case, the magnetic 
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field, through the vector potential A, can be expressed in terms of two scalar functions P and Q 

^ -d+l-^in -Up). +Q} = VP X V0 + QV0 . (8) 



r sin r dO ' 5r 

We remark that the velocity field in cylindrical stationary state accretion-ejection model is a 
three-component field. This is necessary because the jets are generated by magnetocentrifugal 
motion of the planar disk plasma to the axial direction. Here, in this spherical dynamic model, 
the velocity field is a two-component field with vg = 0, because the jets are generated by 
first pulling the disk plasma to the center and then redirecting it to space radially through 
erruptions. 

By self-similar solutions in time, we mean a special class of time-dependent solutions where 
the time and space parts of the physical quantities are in a separable form. The time part is 
described by the evolution function y{t), and the space part will be solved self-consistently by 
separation of variables. The concept of self-similar dynamics is closely related to self-organized 
states that often have minimum energy under given constraints. Having this in mind, we now 
transform the independent variables from {r,9,t) to {r],9,y), and proceed to determine the 
explicit dependence of y in each one of the physical quantities with this radial velocity. First, 
making use of Eq.(7), Eq.(l) renders 

I + ^^.('■'-p) = (| + + "(I + 1) = » • (3") 

Considering the second equahty, the first bracket is the convective time derivative in Euler 
fiuid coordinates, and this amounts to the time derivative in Lagrangian fiuid coordinates. We, 
therefore, have 

dp_dy^ ^dv^2v ^ dp_dy^3pdy ^ ^ ^ q 

dy dt dr r dy dt y dt dy y dt 

Solving this equation for the y dependence gives 

p{f,t) - \p{v,0) . {9b) 

y 



As for Eq.(6), with F = (p/pT) it follows that 
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OF dF dFdy ^ , 



Using the representation of Eq.(8), Eq.(3) is represented by the following two equations 

dP dP dPdy ^ , , 

The first equation gives the function P as 

P{r,t) = \P{V,9) . (116) 

y 

As for the function Q, the right side of the second equation vanishes for rigid rotation to give 

V(j) — a;orsin^ , (13) 



^ ^ f n\ dQdy Qdy dQ Q dy 

Q{f,t) = -^Q{v,e) . (146) 

y 

We note that amounts to a constant plasma rotation, which has to be distinguished from 
the Keplerian disk rotation of the solid material. Comparing to the erruptive time scale of 
our dynamic approach of jets, this rotation rate is negligible and we take cuq — 0. As for the 
magnetic field components, they are given by 

1 - 11 IdP 

Br = -^Briv,e) = +-^-mbbbbbbb-— , 

Bo = \bo{V:0) = , 

y-^ 7] sm a or] 
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1 1 



77 sin 6 
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3. Self-Organization 

Before proceeding with the analysis, let us recall on the fundamentals of self-organiation, 
in particular on MHD systems. We note that MHD equations, like Navier-Stokes equations, 
have quadratic invariants in the absence of dissipations. In MHD systems, there are three of 
them. They are the total energy (plasma and magnetic), magnetic helicity, and cross-helicity. 
Because of the existence of multiple invariants, the system tends to develope self-organized and 
self-similar states through dissipative processes regardless of the details of the initial conditions 
[Hasegawa 1985]. A simple example in fluid mechanics is the developement of shock waves from 
an initial explosion. Another example is that fluid vortex rings (solitons) in air are often formed 
in fast upward drafts of smokes. Should we consider a thin layer of oil heated from below, a grid 
of complex highly organized hexagonal convective cells would develope regardless of the details 
of the initial conditions. We note that, for self-organized states, the memories of the initial 
conditions are lost. In other words, we can not trace a self-organized state backwards in time 
to its initial conditions. They are lost in the dissipative processes that lead to organization. 
The underlying physical arguments for self-similar solutions are also discussed in detail by Low 
[1982a]. For these fundamental reasons, although complex self-similar solutions are only a 
subset of general time-dependent MHD solutions, where most of them are not self-similar, they 
are prone to develope in nature with simple initial configurations. 

Numerically, starting from MHD equations with any fluctuations in a given initial conflg- 
uration, self-organization could be reached since they are insentive to the details in the initial 
conditions. Analytically, the nature of self-similarity implies that the dependent variables (r, t) 
of a physical variable appear in separable form as we have done in Eqs.(9-14). As a conse- 
quence, the solutions will be obtained by the method of separation of variables. Naturally, 
this imposes severe restrictions of the physical system where such a procedure is feasible, such 
as the dimensionality and symmetry. Most of the self-similar solutions arc established in ax- 
isymmetric systems whether cylindrical or spherical. Nevertheless, Gibson and Low [1998] have 
made a great leap in establishing a three-dimensional spherical solution. Under the framework 
of separation of variables, different kinds of solutions can be obtained for the same system de- 
pending on the choice of constants. Likewise in MHD systems, we can have different self-similar 
solutions to account for different phenomena depending on how we separate the constants. In 
the next section, we choose to solve the system with an oscillating radial solution because 
this solution gives magnetic toroids along the jet that match with observations. A monotonic 
decreasing or increasing radial solution in power form is also possible [Lynden-Bell and Boily 
1994]. Although this monotonic solution is not relevant for galactic jets, it could be useful for 
other natural phenomena such as in two-dimensional interplanetary magnetic ropes [Tsui and 
Tavares 2005] and some other astrophysical objects. 
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4. Low Model 

Comparing the equation of F = [p/ p^) and the equation of P, we conclude that F — F{P) 
is a functional of P or F = P(P) ■ We remark that usually we can not make the above statement 
just based on the similarity of the governing equations. It is only possible when we are under 
the framework of self-similarity. We can, therefore, write the 9 dependence in p and p in terms 
of P to get 

P = piv,P) , (15a) 
p = 1-F(rj,e)p'{ri,9) = l-F{P)py(rj, P) = ^7^(77, P) . (156) 

y y y 

Furthermore, the 77 and P dependences should be in a separable form in both 

p(77,P) = PMP2{P) , 

P(V,P) = Pi(v)p2(P) , 

so that, with an adequate 7, F could come out as a functional of P only. Making use of Eq.(4) 
to eliminate the current density in Eq.(2), we get the momentum equation which has three 
components. First, we examine the (f) component which contains only the magnetic force 

di^dQ _ dPdQ ^ ]_ dPdQ _ d^dQ 

Or 89 89 8r ~ y^^ 8r] 89 89 8r]^ ~ ' ^ 

Q = Q{P) . (166) 

The vanishing of the magnetic force in the azimuthal direction described by Eq.(16a) implies 
the functional relationship given by Eq.(16b). As for the 9 component, with the knowledge of 

p = p{r), P)/y^'^, it reads 

dP ^d^P 1 . 8^ 1 8P^ ^8Q ,8P 2.2n9PidP^ 
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We remark that the first three terms of this equation represent the force-free field equation for 
JxB = [My 1984, Low 1986, Low and Lou 1990, Lynden-Bell and Boily 1994]. The last 
term is the plasma pressure. This equation would be independent of the evolution function y 
should we consider the 7 = 4/3 case pioneered by Low. 

We note that a nonhnear equation like Eq.(17a) has to be solved subject to the boundary 
conditions of a given physical problem. For example, should the problem on hand has an in- 
finite external domain and power form radial solutions are physically reasonable, then we use 
fractional powers of P for the functionals of Q and p [Lynden-Bell and Boily 1994]. In our 
present case, we are interested in decaying oscillating solutions in 77. We, therefore, write 

Q{P) = aP , (18a) 



P{V:P) = pMP2{P) = iv-W'P' + C) , (186) 

where (7 is a positive constant independent of coordinates and the functional P. This is the 
simplest representation where the first equation gives a linear dependence of P and the second 
equation reflects the positive definite nature of plasma pressure. This choice of constant C 
in plasma pressure profile differs from the spherically symmetric radial coordinate dependent 
additive term in Eqs.(21) and (22) of Low [1984b]. This difference of representation stems 
from the view that Low considers Eq.(17a) as an equation that solves for the plasma pressure 
under a given P. For this reason, the homogeneous solution corresponds to the spherically 
symmetric gasodynamic solution. We regard Eq.(17a) as an equation that solves for P under 
a given plasma pressure that has a positive definite separable form as in Eq.(18b). The choice 
of Eq.(18a) gives dQ/dP = a, a constant. With 6^ = 2//6'^, Eq.(17a) now reads 

^'^ + sm9^{-^^)+rj'a'P + b'sin'eP = . (176) 

Writing P{ri,9) = R{r])Q{9), x — cos9, and with n{n + 1) as separation constant, Eq.(17b) 
becomes 

rf^ + [a^rf -n{n+l)]R = , 
(l-x2)^^^^ + [n(n + l) + 62(l-x2)]e(x) =0 . 
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The R{r)) equation can be solved readily to give 

R{V) = {avY^^Jn+i/2{ar)) . (19a) 

Such a spherical Bessel functional was used by Low [1984b] in Eq.(28) of his paper as one of 
the numerical examples in spherical two-dimensional self-similar MHD to model coronal mass 
ejections. 

The Q{x) equation with finite plasma pressure, ^ 0, can be solved by power series 

e{x) = J2 ' (19^) 

(m + 2){m + 1)0^+2 = [m(m - 1) - n{n + 1) - &^]am + , 
603 = -[n{n + 1) + b'^]ai , 

2a2 = -[n{n + 1) + &^]ao , 

where the infinite sum in Eq.(19b) starts from m — 0. There are two independent solutions. The 
first one has oq = ©(0) 7^ and oi = 0, and the second one has ao — and oi = dQ{0)/dx ^ 0. 
They correspond to even and odd powers of the series respectively 
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In the absence of plasma pressure with 6^ = 0, the above solution reduces to 

e{x) = (l-x2)^^l^ = -n{n + l) / Pn{x)dx , (196') 
dx J I 

where Pn{x) is the Legendre polynomial. The even solution corresponds to n odd, and the odd 
solution is otherwise. 
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5. Self-Similar Magnetic Field 

With the solution P{r],9) — R{r])Q{9) established by Eqs.(19), the magnetic field compo- 
nents are 

77"^ ax 
1 ldR{r]) 1 

We note that the self-similar evolution of the MHD plasma distorts the dipole-hke magnetic 
field by generating an azimuthal component of the magnetic field with finite a^. To grasp the 
magnetic structure given by Eqs.(19a) and (19b), we first examine the case with 6^ = 0. In 
this special case, the even and odd power series of Eq.(19b) will terminate at a finite number 
of terms when m = n -|- 1 to give Eq.(19b') and Q{x) = at the poles x — ±1. 

The self-similar radial structure R{r]) given by Eq.(19a) allows oscillations in t] if at] is 
sufficiently large, which means that the azimuthal magnetic field is sufficiently large. The 
meridian structure 0(x) given by Eq.(19b') also oscillates in x. Let us denote rji and xj as 
where R{i]) and Q{x) vanish. We remark that rji are circles of constant 1], and xj are spokes of 
constant x = cos 6. Consequently, {r]i,Xj) divide the [rj — x) plane in many smaller regions. On 
r/j, we have 5^ = and = 0, with B,. and B^ changing signs across rji. The only component 
that does not vanish completely on r/i is Bq. Referring to the complete expression of Bq above, 
this magnetic component is modulated by Q{x) so that it changes sign on the circle of constant 
1] on crossing each spoke region of Xj. On xj, we have Bg = and B^ = 0, with Bg and B^ 
changing signs across Xj. The only component that does not vanish completely on xj is Br- 
Referring to the complete expression of Br above, this magnetic component is modulated by 
R{i]) so that it changes sign on the spoke of constant x on crossing each circular region of rii. As a 
result, axisymmetric closed magnetic field lines are formed in these regions generating toroidal 
belt plasmoids circumscribing the z-axis of symmetry. We call these secondary plasmoids, 
and call the larger self-similar plasmoid embedding all the secondary plasmoids the primary 
plasmoid. Neighboring plasmoids have field lines circling in opposite sense. If one plasmoid 
has field lines in clockwise direction, the adjacent one has them in counter clockwise direction. 
The azimuthal components also rotate against each other. 

In each region bounded by {r]i,r]i+i) and {xj,Xj+i), the topological center defined by 
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dR{ri)/dri = and dQ{x)/dx = has Br = and Bg = 0. This is the magnetic axis of 
each toroid. The field fines about tfiis center are given by 

f: = ^ = . (20«) 

dr rdO r sm Odq) 

By axisymmetry, tfie magnetic field components are independent of 0. For this reason, the third 
group of the above equation is decoupled from the first two groups. The B^ field circles about 
the z-axis of symmetry without twisting. It is simply superimposed on the B^ — Bg field lines. In 
terms of Fourier components e^™''^', this means m = 0. For the field lines on an (r — 9) plane, we 
consider the first equality between Bj. and Bg which gives P — R{r])Q{x) equals to a constant or 

P{r),x) = {ar))'/^Jn+i/2{avMx) = C . (206) 



In other words, the nested field fines are given by the contours of P{ri, x) on the (r — 9) plane. 

As an example, the field lines for n = 4 and 6^ = arc shown in Fig.(l) for axisymmetric 
secondary plasmoids. We have taken ariQ = 8.2 at the first zero of the spherical Bessel function 
where rjo is the radial label of the plasma boundary. The horizontal axis isa; = 0or^^ = 7r/2 
and the vertical axis is a; = +1 or = 0. For n = 4, Q{x) vanishes at x = 0, ±(3/7)-'^/^, ±1, 
and R{i]) vanishes at arj = 0.0, 8.2. These arc locations where Bg = and Br = respectively 
dividing the quadrant in two regions, as indicated in Fig.(l). These dividing lines arc obtained 
by solving the contours with C = numerically. The radial grids are equally spaced in 1]. The 
meridian grids are equally spaced in x. By converting to 6 through x = cosO, it generates 
an uneven grid distribution in 6 that appears in Fig.(l). Closed field lines are also shown in 
each region. Negative value contours of C = —0.3, —0.5, —0.7 show the field lines in the region 
adjacent to the axis of a; = 0, and positive value contours of C = +0.3, +0.5 show the field lines 
in the region adjacent to the axis of a; = +1. At the topological center of each region, we have 
Q{x) maximum and R{r)) maximum, thereby giving Bg = and i?,. = at the same location. 
Since 



27rr sm9B^ = 2TraP , 



the center has the maximum of line integral of B^ about the axis of symmetry. Adjacent 
secondary plasmoids have oposite signs of i?^. Should we take arjo = 11.7 as the second zero 
of the spherical Bessel function, one additional shell of secondary plasmoids would be added 
to Fig.(l). Furthermore, because of the (1 — x"^) factor in Eq.(19b'), Q{x) is null at the poles 
with X — ±1 where the magnetic field is purely radial. If n ^ 1, the null at the poles is tightly 
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surrounded by a polar lobe. Along the polar lobe, the embedded belt plasmoids are tightly 
wounded about the polar axis, and we name them secondary plasmoids. 

We are now ready to examine the general case of 6^ 7^ 0. This case of 6^ 7^ is known as 
the finite /3 case where (3 is the ratio of plasma pressure to magnetic pressure. The presence of 
the plasma pressure with finite 6^ distorts the linear solution of Q{x). With finite /3, the series 
usually has a finite value at x = ±1 which leads to singular magnetic fields there because of 
the (1 — x^y/'^ factor in the denominator. Constrained by nature to regular solutions, 6^ has to 
be the eigenvalues such that Q{x) remains null &X, x — ±1. To search for these eigenvalues, we 
evaluate ©(x), Eq.(19b), at the boundary x — for different values of 6^ with a given n. The 
eigenvalues are shown in Fig.(2) as the intercepts of ©(+1) = at 119, 261, 425 for even series 
with n = 13. As for n = 23, they are 200, 422. The even series has 0(0) 7^ and de{0)/dx = 
and the odd series has ©(0) = and dQ{0)/dx ^ 0. The even eigenfunctions with n = 13 and 
= 119, and with n = 23 and }? = 200, are shown in Fig.(3). With n = 13 and 6^ = 119, 
there are seven nodes in the inverval (0, 1). Should we take the second eigenvalue = 261, one 
more node would be added. With n = 23 and 6^ = 200, there are twelve node in (0, 1). The 
next eigenvalue — 422 would add one more node as well. 

The fact that the plasma pressure has to be as such that it is the eigenvalue of the Q{x) 
equation appears to be a very restrictive constraint for the self-similar solutions. Nevertheless, 
we note that the plasma pressure appears in the ©(x) equation where the separation constant 
n{n + 1) is an as yet unspecified free parameter. As a result, there will be an adequate n for 
almost any given plasma pressure. For example, with &^ = 422, instead of being the second 
eigenvalue oi n — 23, it could be the first eigenvalue of some different n that is larger than 23. 
The important point is that for any plasma pressure b^, it will coincide to one of the eigenvalues 
of some n. The chances are that n is a large integer which, as we will show in the next section, 
leads to good jet coUimation. 
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6. Jet CoUimation and Vortex Structure 

One important point we should point out is that, although we have taken a spherically 
symmetric radial expansion in the plasma velocity, the plasma density and pressure need not 
be symmetric, and much less the magnetic fields. The plasma profiles will be discussed later. 
Here, we discuss the magnetic fields that are given in terms of R{ri) and Q{x) as hsted in the 
preceeding section. The maximum of Q{x) at x — gives Br — on the equatorial acreation 
disk plasma where the magnetic field has only a Bg meridian and a B^ azimuthal component. 
Off the disk, due to the nodes in Fig. (3), there are also maxima in space, so that the magnetic 
fields are structured in lobes other than the equatorial lobe. 

To consider the collimation of polar jets, we examine the magnetic field components. For 
Bg and B^ components, we plot the function 0(x)/(l — x^)^/^ in Fig. (4) with the corresponding 
eigenvalue 6^ = 119 for n = 13 according to Fig. (2). It peaks up off the scale to (-3.0) as x 
approaches unity, but plunges to null at x = +1. As n increases, the peak edges closer to 
X = +1. This is clearly shown in Fig. (4) with = 200 for n = 23. The peak then goes 
to (+3.8) in this case. As a consequence, the peak rises up in amplitude but narrows down 
in width. To examine Br, we plot the function dQ{x)/dx in Fig. (5) which shows that the 
magnetic field is purely radial on the polar axis. As we have said, this radial magnetic field 
changes sign as each region of rii is crossed. With n = 13 and n = 23, the peak on the polar 
axis is (+71) and (-155) respectively. The magnetic energy is the quadratic quantity of this, 
and it peaks even more with respect to the off axis lobes. The detail structures near the polar 
axis are shown in Fig. (6) for Bg and B^j,, and in Fig. (7) for Br- The magnetic field, therefore, 
converges to the polar axis as n increases leading to a jet structure. Comparing the equatorial 
lobe of the magnetic field with the polar lobe, it is clear that the polar lobe is much narrower 
than the equatorial one. This lobe pattern is similar to directional antenna arrays. Due to 
the oscillating nature of Eq.(19a), there are spherical toroids, or secondary plasmoids, formed 
by closed magnetic field lines along this narrow peak. The number of secondary plasmoids 
embedded in this peak measures the number of zeros of Eq.(19a) contained within arjQ. The 
magnitude of arjQ also indicates the large amplitude of B^ through Q = aP of Eq.(18a) which 
is itself consistent to the magnetic field convergence onto the polar axis. 

To understand the magnetic structure of the jet, we plot the contours of P given by 
Eq.(20b) which amount to closed field fines of embedded secondary plasmoids. A contour plot 
in Fig. (8), similar to Fig.(l) but with n = 13 and — 119, shows plasmoids along the polar 
cone that goes from x — 0.96 to x — 1, as in Fig. (7), which corresponds to a cone angle of about 
14 degrees. The radial label goes from arj — to arj — 40, and this interval is divided into 
seven regions set by the zeros of the spherical Bessel function. The first zero is at arj — 18.40 
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and the subsequent zeros are approximately equally spaced, as indicated in Fig. (8). Negative 
contours of C = —0.010 are plotted in the first, third, fifth, and seventh regions only. Likewise, 
positive contours of C = +0.005 are plotted in second, fourth, and sixth regions. To avoid over 
crowdedness, the C = —0.005 contours are not plotted in the odd numbered regions, and the 
C = +0.010 ones are not shown in the even numbered regions. Once more, we recall that the 
magnetic field lines of these secondary plasmoids circle in alternating sequence. The boundary 
of the jet is given by the root of Q{x) with xj = 0.96. At this boundary, the magnetic field is 
purely radial but alternates in sign. The other boundary is at xj = 1.0 where the magnetic field 
is also radial and alternates in sign opposite to that of Xj = 0.96. Together with the roots of 
R{i]) at r]i where the magnetic field is meridian, the plasmoids are magnetic vortices bounded 
by closed field lines at the border of each region. By axisymmetry, these poloidal magnetic 
contours rotate about the polar axis to generate tightly wounded toroids. Furthermore by 
axisymmetry, the toroidal azimuthal field lines are decoupled from the bipolar poloidal field 
lines. The field lines are, therefore, two-dimensional. As in the discussions concerning Fig.(l), 
this amounts to the m = mode in the Fourier transform e*'"'^ of the toroidal dependence. 

The magnetic structure of this spherical temporal self-similar model, which shows a se- 
quence of magnetic toroids along the jet, differs from the continuous helical field line structure 
in the cylindrical spatial self-similar model, complemented by numerical simulations [Casse 
2004]. To understand the differences, the spatial cylindrical model has an imposed initial mag- 
netic configuration everywhere in space. Jets are formed by transporting the disk plasma from 
the disk plane to the axis through magnetocentrifugal action. They are maintained along the 
axis by bringing disk plasma continuously in steady state to overcome transport losses. In this 
description, spatially self-similar jets in z/r were formed in their present position in the distant 
past, and are maintained there continuously at the present and in the future. 

As for the temporal self-similar spherical model, the disk plasma is first accreted to the 
galactic nucleus which builds up the plasma pressure there. While the plasma is in this bounded 
region, self-organization is nourished through dissipations to self-similar structures in 77 = 
r{t)/y{t) calculated here, including the jets. The force-free configuration and the vortex nature 
of the magnetic field are akin to the quadratic quantities of magnetic helicity and cross helicity 
of the ideal MHD system. We note that rj is time invariant, such that the temporal self- 
similar 77 profiles in magnetic field and plasma parameters maintain their forms at the very 
beginning. Later on, the self-similar MHD plasma errupts from the center outward, and this 
entire structure expands into the previously void space as time progresses, and continues to do 
so until the plasma density and magnetic field intensity get so low that they fade into space 
light years away. The bounded state and the erruption are described by the time evolution 
function y{t), and it will be solved consistently in the next section. 
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7. Evolution Function 



The r component of the momentum equation reads 



,dv dv, d , GM 



' dt dr dr ^ ' 

/X r sm fc* or or^ o9 sm 6 oO oP 

The term dp/dr on the left side refers to radial derivatives for the explicit dependence and 
the implicit dependence in P. Making use of the 9 component, Eq.(15a), the equation above 
becomes 

.dv dt' d , GM 
p\-T- + v—\ + — p(r, P(r, e)) + p- 



2 



dt dr dr r 



^2^^r 2 • 2/)^Pi 9P dp 



/J, rsinO dr dP dr dP 

The right side is just the radial derivative of plasma pressure on the implicit dependence in 
P. This cancels the corresponding term on the left side leaving only the explicit radial derivative 

dv ^dv dp{r, P) GM _ ^ 
^ dt ^ dr dr ^ 

In terms of self-similar parameters, this equation reads 

oA l^.^ldp GM^ .d^y lAdp GM ^ , , 

dt^ r) pdrj rj'' dt^ rj pdrj rj'' 

With a as the separation constant, and H as an integration constant, the evolution function 
is, therefore, described by 

l-2(ir+y]=H. (23a) 
| = ±[2(ff-^)]"V (23(,) 
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To understand the meaning of a, we note that plasma acceleration in Lagrangian coordi- 
nate is 

dv dPy ar) arf 
dt ^ dt"^ 

A negative a means an outward decelerating flow, or an inward accelerating flow. The deceler- 
ation gets smaller as y or as r gets larger, the acceleration gets larger as y or as r gets smaller. 
As for the meaning of if, we take the hmit y — oo which gives H/At^ — [v'^ /2)/ Ai^rf by using 
r — rjy. Physically, this is the asymptotic radial kinetic energy per unit mass per unit area of 
a Lagrangian fluid element. For circular orbits, we would have H — 0. Furthermore, we can 
rewrite Eq.(23a) making use the above physical expression of a to get 

prj^H = ^pv^ + (r • F)r , 

where F denotes all the forces on the right side of Eq.(2). It is clear that H measures the total 
energy of the fluid element. 

To consider the accretion disk magnetospheric plasma, we take a — — |a| < negative for 
an inward axcelerating flow. With H — as the boundary condition at infinity, we get 

dt S y 

The two signs for the square root correspond to the outward and inward flows. Taking the 
negative sign for the inward flow, the equatorial accretion disk plasma starts at infinity with 
an asymptotically zero radial velocity and ends near the center with a large influx. 

As to describe the polar jets, we again take a = —\a\ < negative. We consider H = ±|if| 
which lead to 

^ = ±[2(±|i/| + M)]i/2 . (246) 
dt y 

For H = +\H\ > and taking the positive sign on the square root, this gives a decelerated 
outward flow. The deceleration gets smaller as y or as r gets larger, with the terminal velocity 
given hy V = -qdy/dt = 77 (2i/)^/^. For H = —\H\ < 0, the outward flow would stop at 
y — \a\/\H\ where dy/dt — 0. This H — —\H\ < outward flow would be followed by an 



(24a) 
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accelerated inward flow, should we take the negative sign on the square root, so that the system 
would be periodic and bounded. We, therefore, see that the system would make a transition 
from a bounded to an unbounded state when H goes from negative to positive. The bounded 
state allows us to define the boundaries of the radial label (0,r7o) such that (0 < 77 < rjo). The 
solutions of the evolution function in Eq.(24b) suggest that jets are a result of an erruptive 
process. This process is fed by plasma accretion from the galactic disk. As pressure builds up 
at the center, the plasma begins to oscillate, or to pulse, as described by the bounded solution 
with H = —\H\ < 0. We believe that it is during this phase that self-organized and self-similar 
structures in plasma parameters and magnetic fields would be nourished through dissipations. 
The dynamics is given by self-similar solutions where the magnetic fields are self-consistent to 
the plasma pressure that plays the role of the eigenvalue. The jet structure emerges along the 
polar axis as the eigenfunction of the magnetic field. By further pressure built-up and energy 
infiux, H = +\H\ > becomes positive, and the structure eventually errupts. 
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8. Rayleigh-Benard Cells 

Now we have presented the general structure of galactic jets. It begins with an influx of 
plasmas from the accretion disk which drives up the plasma pressure at the galactic nucleus. 
The MHD plasma responses to self-organization by oscillating, or pulsating, periodically in 
time. As H goes from negative to positive, the periodic mode goes to an erruptive mode. 
The result of self-organization is a structured configuration in space. This structure is not 
unique. For galactic jets, we have solved the system with the method of separation of variables 
as such that it ressembles to observations. This structure contains basically convective cells 
in space where magnetic field lines in adjacent cells rotate in opposite sense. Although such 
a structured configuration with consistent evolution function is in accordance with the time- 
dependent MHD equations, there is no mention of the initial configuration that could lead 
to such self-similar states. Consequently, such a highly organized complex structure could be 
regarded as an artificial result due to special mathematical constructions, which might not have 
any relevance to the physical jet system. In order to bring more reality to this analytic result, 
we recall the classical case of Rayleigh-Benard fluid self-organization. Consider a thin layer 
of oil heated from below, observations tell us that this simple homogeneous configuration will 
develope an array of identical hexagonal convective cells if the temperature gradient across the 
layer is sufficiently large. The velocity streamlines of adjacent cells rotate in opposite sense. 
Should the initial state of this oil layer be altered by arbitrary fluctuations, same convective cells 
would still appear after the initial fluctuations are dissipated. Numerically, such self-organized 
complex structure could be reached from the simple homogeneous initial configuration if the 
code is adequately pushed in the correct direction. 

Guided by this fluid example, the self-similar jet structure could be regarded as the 
Rayleigh-Bernard equivalent in the MHD system. We could choose an initial state by tak- 
ing n = 1 in our solutions with plasma bounded between an outer sphere rjQ and an inner 
sphere rj^,. This would give a large global long wavelength structure in the spherical layer. We 
consider an equilibrium state with a = —\a\ < and H = —\H\ < with y = \a\/\H\ and 
dy/dt = 0. Let us take this moment as t = 0. An energy source in terms of pressure is supplied 
by the accretion disk to pump the MHD plasma at the lower boundary rj^,. At some moment, 
y begins to depart from its equilibrium to fall towards the center due to perturbations. Like 
in the oil layer case, this could trigger convections so that the long wavelength global struc- 
ture cascades to short wavelength structures accepted by the system. Cenvective cell scaling 
becomes smaller as n gets larger. Different from the flat layer fluid case, the convective cells in 
this spherical MHD layer are not identical, since there is a focusing effect to the polar axis. 
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9. Mass Density Profile 

As for the spatial part, we have 

Idp GM 
p orj T}'^ 

With p given by Eq.(18b), the above equation gives 

= {r'){-^W'P' + C) = pMP2{P) . (25) 

In our present case, we are considering a spherical shell domain that the radial label r] = is 
excluded. The negative powers of rj in plasma pressure and mass density do not cause singu- 
larity. We have set a = 0, with vanishing net force on the flow, to obtain the second equality 
in Eq.(25) to be compatible with the functional form of F = F{P) with 7 = 4/3. Nevertheless, 
we could relax this condition to a ~ such that ar]^ « GM. To understand the implication 
of such inequality, we remark that 

^ _ dv/dt 

GM ~ GM/t^ ^ ■ ^^^^ 

Therefore, it is just the acceleration of the flow to that of the central mass gravitational acceler- 
ation. Since dv/dt gets smaller as r gets larger with 1/r^ scahng, this condition estabhshes that 
self-similar configurations can be organized, and bounded oscillations of the H — —\H\ < 
case can take place, as long as the plasma acceleration is much less that the gravitational ac- 
celeration. To close the entire self-similar system, we now come to the functional F{P) for the 
adiabatic equation of state. With the results in Eq.(18b) and Eq.(25) using a 0, we have 

m = = i^f'p-^"'ip) = (^f'(b^p' + . (27) 



To understand the plasma structure of the jet, we note that the plasma pressure and mass 
density are b'^P'^{ri,x) + C dependent which means positive definite + C dependent, ac- 

cording to Eq.(18b) and Eq.(25) respectively. The plasma pressure and the mass density have 
their minimum at P{r),x) — {ar)y^'^Jn+i/2{0'V)^{^) = 0- This minimum is positive nonzero 
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because of the positive integration constant C in Eq.(18b). To discuss the plasma structure in 
the jets, we reproduce the part of Fig. (3) in the range 0.94 < x < 1.0 in Fig. (9). The segmented 
hne for n = 13 has 0(x) = at the cone boundary x = 0.96 and at the cone center x = 1.0. 
In between, there is a minimum at a; = 0.985. Because of the quadratic dependence G^(a;), 
these features correspond to a cavity-hke structure for the plasma pressure and mass density of 
the jet, compatible to numerical simulations of the time-dependent dissipative MHD equations 
[Casse 2004, Zanni et. al. 2004]. Other than the zeros of Q{x), P{f],x) also vanishes at radial 
locations where ^,1+1/2(0^7) = 0. With the quadratic dependence 

there are ripples for plasma pressure and mass density along the radial label. The radial 
function has two contributions. The part of the Bessel function Jn+i/2{(^Vl) has an oscillating 
nature with decreasing amplitude, while the part of {arj) has an increasing amplitude that 
helps to maintain the ripple level of R^{r])- The peaks of these ripples are at the topological 
center of each plasmoid which is the magnetic axis of the toroid where has a maximum 
line integral. These ripples are also seen in numerical simulations [Zanni ct. al. 2004]. It is 
quite surprising that the results of this spherical temporal self-similar MHD model agree rather 
well qualitatively with numerical simulations of the time-dependent MHD equations for the 
cylindrical magnetoccntrifugaling model. Both of them give cavity structures in the transverse 
direction and ripple structures in the longitudinal direction. 

The polar jets, therefore, have a periodic and approximately equally spaced concentration 
of plasmoids in its longitudinal direction, except the first region which is more extensive. The 
closed magnetic field lines for adjacent plasmoids rotate in opposite sense. The plasma pressure 
and mass density have a hollow conic structure with ripples along the radial direction. Such 
periodic structure happens to be compatible with radio observations where high intensity islands 
are mapped. These islands are usually thought to be periodic ejections of mass from the 
accretion disk. According to our model, they are rather the internal spatial arrangments of an 
expanding jet driven by an erruptive event. 
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10. Discussions and Conclusions 

One of the main objections of self-organized plasmoid representation in free space in the 
absence of an adequate boundary is that it apparently violates the Virial theorem which states 
that [Schmidt 1966] 

l^ + J ^^^^^ ^2{T + U) + W'= + W''- J XkiPik + T,k)dS , (28) 

where / is the moment of inertia of the plasmoid, G is the momentum density of the electro- 
magnetic field, T and U are the kinetic and thermal energies of the plasma, and W'^'^ are 
the electric and magnetic energies in the volume, Pik and Tik are the plasma and electromag- 
netic stress tensors. Taking the volume to cover the entire plasma and field, the surface term 
on the right side vanishes. In laboratory plasmas, this surface could be the machine vessel. 
Should the plasmoid be stationary, the volume term on the left side would be null, and the mo- 
ment of inertia would be accelerating since the terms on the right side are all positive definite. 
While this statement has no confiict with the unbounded solutions, it apparently contradicts 
the stationary state of the bounded solutions. Nevertheless, this argument has overlooked the 
asymptotically bounded nature of the H = —\H\ < plasmoid state. In this asymptotic case, 
we have dl /dt = {dl / dy){dy / dt) = so that / is stationary. The fact that d^I /dt^ > implies 
that / is at an asymptotic minimum, not an acceleration of /, which complies with the Virial 
theorem. 

The classical accretion-ejection model of Blandford and Payne [1982] is a time-independent 
stationary state model with spatial self-similar MHD solutions in cyhndrical geometry (r, 0, z) 
with Alfvenic plasma flow velocity plus a rotation, all with Keplerian scaling. In this model, the 
jets are formed by convecting the magnetospheric disk plasmas from the disk plane to the axis 
by magnetocentrifugal action through the magnetic field lines with low inclination angles to 
the disk plane. The angular momentum of the plasma on the disk plane is focused to the axis. 
The jets were put in place in the distant past according to the spatial self-similar solutions in 
z/r, and are maintained there by continuously transporting disk plasmas to the axis to sustain 
the axial outflow. CoUimation to the axis is accomplished by magnetic hoop force. 

Here, we have taken a dynamic view where jets are the consequences of erruptive events, 
based on time-dependent MHD equations in spherical geometry {r,9,(f)). The disk plasmas 
are accreted to the galactic nucleus where plasma pressure is built up to cause an erruption. 
The radially symmetric expanding velocity interacts the plasma with the magnetic field self- 
consistently through the MHD equations to generate spatial structures. Due to the existence of 
multiple quadratic invariants in the absence of dissipations, MHD systems have the tendency 
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of developing self-organized and self-similar states through dissipative processes. The force-free 
configuration and the vortex nature of the magnetic field are akin to the quadratic invariants 
of the magnetic helicity and cross helicity of the MHD system. For these reasons, although 
temporally self-similar solutions are only a subset of general time-dependent MHD solutions, 
these self-similar configurations are prone to develope in natural phenomena. 

We, therefore, describe these spatial structures by self-similar temporal solutions in 77 = 
r/y. In this self-similar spherical model, consistent self-similar representations of plasma pres- 
sure, mass density, and magnetic fields are worked out for the 7 = 4/3 BC Low model, with 
special emphasis on the finite (3 case. The spatial distribution of the magnetic field is described 
by an equation where the plasma pressure acts as the eigenvalue. The nature of this eigenvalue 
equation is as such that the magnetic field gets converged to the polar axis as a response to the 
high plasma pressure directly related to the eigenvalue. Since the separation constant n{n-\-l) 
in this eigenvalue equation is a free parameter that as yet to be specified, there will be an ade- 
quate n for almost any given plasma pressure. Although the radial plasma velocity is isotropic, 
the spatial structures are not. They could be highly coUimated along the axial direction, and 
expand into the previously void space as time progresses. 

Although other types of solutions are permitted, we have specifically examined the solutions 
that bear resemblance with jet features with n ^ 1 and 0770 ^ 1- Since plasma and magnetic 
field are frozen into each other, plasma outflow is also coUimated to the polar axis. The 
existence of regions of closed fleld lines along the primary polar magnetic lobe permits secondary 
plasmoids be embedded in it. These secondary plasmoids appear to be compatible to the 
observed islands of radio intensities. The time evolution function of the radial velocity consistent 
to the temporal self-similar solutions has different types of solutions according to the sign of 
H. Although the equatorial disk magnetospheric plasma is not addressed here, the accelerated 
accretion of plasma influx could be modelled with if = as the boundary condition at inflnity. 
As for the polar jets, H < gives a bounded oscillating, or pulsating, solution. This bounded 
stage, we beheve, nourishes the self-organized and self-similar states. With if > 0, it gives 
an unbounded expanding solution with a high terminal velocity. It is apparent that plasma 
pressure due to accretion is the prime driving force that determines the value of H. The 
bounded oscillation would make a transition to the unbounded expansion as H goes from 
negative to positive. According to our model, jet structures are, therefore, considered as a 
spatial conflguration that has been expanding continuously into space. 
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Fig. 1. — The field lines on the (r — 6) plane are plotted in one quadrant with arjo = 8.2 as the 
first zero of the spherical Bessel function. The quadrant is divided into two regions by C = 
with plasmoids having C = —0.3, —0.5, —0.7 and C = +0.3, +0.5. 
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Fig. 2. — The function Q{x) at x = +1 is plotted against the parameter 6^ to identify the 
eigenvalues for n = 13 in segmented line and n = 23 in solid line. 




Fig. 3. — The eigenfunction Q{x) is plotted against x with eigenvalue — 119 for n = 13 in 
segmented line and — 200 for n = 23 in solid line to show its dependence on n. 
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Fig. 4. — The function B(x)/(1 — x^)^/^ is plotted against x with eigenvalue 6^ = 119 for n = 13 
in segmented line and = 200 for n = 23 in solid line to show the structure of Bg and in 
space. 




Fig. 5. — The function dQ{x)/dx is plotted against x with eigenvalue 6^ = 119 for n = 13 in 
segmented line and = 200 for = 23 in solid line to show the structure of Br in space. 
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Fig. 6. — The function 0(x)/(l — x^)^/^ is plotted against 0.94 < x < 1.0 with eigenvalue 
6^ = 119 for n = 13 in segmented line and 6^ = 200 for n = 23 in solid line to show the detail 
structure of Bg and near the polar axis. 




Fig. 7. — The function dQ{x)/dx is plotted against 0.94 < x < 1.0 with eigenvalue 6^ = 119 
for 72 = 13 in segmented line and 6^ = 200 for n = 23 in solid line to show the detail structure 
of Bj. near the polar axis. 
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Fig. 8. — The field lines on the (r — 9) cone around the polar axis are plotted with arjQ = 40 and 
a not so large n — 13 for better resolution. The regions are divided by C = with plasmoids 
having C = -0.010 and C = +0.005. 
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Fig. 9. — The eigenfunction G(x) is plotted against 0.94 < z < 1.0 with eigenvalue IP' = 119 
for n = 13 in segmented line and 6^ = 200 for n = 23 in solid line to show its dependence on n. 



